{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 傅里叶变换  \n",
    "> 仅帮助自己理解使用，并非严谨的数学证明。摘自知乎。  \n",
    "\n",
    "$\\int^{t_0+T}_{t_0}\\cos{x}\\,dx\\,=\\,-\\sin{x}|^{t_0+T}_{t_0}$ 因为$\\sin{x}$是周期为$T$的函数，所以积分的值为$0$。  \n",
    "\n",
    "### 三角函数的正交性  \n",
    "对于一个`三角函数系`：$1,\\,\\cos{x},\\,\\sin{x},\\,\\cos{(2x)},\\,\\sin{(2x)},\\,\\ldots,\\,\\cos{(nx)},\\,\\sin{(nx)}$，其中任何两个不同函数的积在区间$[t_0, t_0+T]$内的积分为$0$。$T=2\\pi$。  \n",
    "\n",
    "#### 三角级数  \n",
    "给定一个初始角频率$\\omega_0$对应周期为$T_0$，之后的角频率必须为$\\omega_0$的正整数倍，相应的周期也应缩短为$\\frac{T_0}{n}$，当$n\\rightarrow+\\infty$时，对应的周期趋近于$0$。\n",
    "\n",
    "### 傅里叶级数  \n",
    "对于一个周期为$T$的函数$f(t)$有：  \n",
    "\n",
    "$$\\left.\\begin{array}{lcc}\n",
    "    f(t) = & \\underbrace{A_0} & +\\,\\sum\\limits_{n=1}^{\\infty} {\\underbrace{[a_n\\cos{(n\\omega t)}+b_0\\sin{(n\\omega t)]}}}  \\\\ \\\\\n",
    "    & 直流分量 & \\cos{(n\\omega t + \\phi_n)}\n",
    "\\end{array}\\right\\}相当于用无穷多个正弦函数逼近$$\n",
    "\n",
    "\n",
    "1. 对于直流分量。$\\int_{t_0}^{t_0+T}{f(t)}\\,dt = \\int_{t_0}^{t_0+T}{A_0}\\,dt = A_0T$，故$A_0=\\frac{\\int_{t_0}^{t_0+T}{f(t)}\\,dt}{T}$。令$a_0=2A_0$，则$a_0=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)}\\,dt$    \n",
    "\n",
    "2. 对于交流分量的第一部分$\\cos$。$int_{t_0}^{t_0+T}{f(t)\\cos{(n\\omega t)}}\\,dt = \\int_{t_0}^{t_0+T}{a_n\\cos^2{(n\\omega t)}}\\,dt=\\frac{a_n}{2}\\int_{t_0}^{t_0+T}{[1+\\underline{\\cos{(2n\\omega t)}}]}\\,dt=\\frac{a_n}{2}T$，  \n",
    "  故$a_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\cos{(n\\omega t)}}\\,dt$    \n",
    "  \n",
    "3. 同理：$b_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\sin{(n\\omega t)}}\\,dt$  \n",
    "\n",
    "可得： \n",
    "$$\\begin{array}{lccccccc}\n",
    "    f(t) = & \\underbrace{\\frac{a_0}{2}} & + & \\sum\\limits_{n=1}^{\\infty}[ & \\underbrace{a_n\\cos{(n\\omega t)}} & + & \\underbrace{b_0\\sin{(n\\omega t)}}&]  \\\\ \\\\  \n",
    "    & a_0=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)}\\,dt &&& a_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\cos{(n\\omega t)}}\\,dt && b_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\sin{(n\\omega t)}}\\,dt &\n",
    "\\end{array}\\tag{1}$$\n",
    "\n",
    "于是，对于周期函数$f(t)$有：\n",
    "$$\\begin{array}{ll}\n",
    "    f(t) = \\frac{a_0}{2} + \\sum\\limits_{n=1}^{\\infty}[a_n\\cos{(n\\omega t)} + b_0\\sin{(n\\omega t)}]  &  其中：\n",
    "    \\left\\{ \\begin{array}{l}\n",
    "        a_0=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)}\\,dt \\\\ \\\\\n",
    "        a_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\cos{(n\\omega t)}}\\,dt \\\\ \\\\\n",
    "        b_n=\\frac{2}{T}\\int_{t_0}^{t_0+T}{f(t)\\sin{(n\\omega t)}}\\,dt\n",
    "    \\end{array} \\right.\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 欧拉公式  \n",
    "> 通过欧拉公式可以把累加写成积分的形式  \n",
    "\n",
    "欧拉公式的定义如下：  \n",
    "$$e^{i\\theta}=\\cos{\\theta}+i\\sin{\\theta}\\, \n",
    "\\left\\{\\begin{array}{l}\n",
    "    \\cos{\\theta} = \\frac{e^{i\\theta}+e^{-i\\theta}}{2} \\\\ \\\\ \n",
    "    \\sin{\\theta} = \\frac{e^{i\\theta}-e^{-i\\theta}}{2i} = -i\\frac{e^{i\\theta}+e^{-i\\theta}}{2}\n",
    "\\end{array}\\right.$$  \n",
    "用欧拉公式替换公式$(1)$中的三角函数，有：  \n",
    "$$\\begin{array}{ll}\n",
    "f(t) & = \\frac{a_0}{2} + \\sum\\limits_{n=1}^{\\infty}[a_n \\frac{e^{in\\omega t}+e^{-in\\omega t}}{2} + (-i)b_n \\frac{e^{in\\omega t}-e^{-in\\omega t}}{2}], \\, 变形 \\\\ \\\\\n",
    "     & = \\frac{a_0}{2} + \\sum\\limits_{n=1}^{\\infty}[ \\frac{a_n - b_n}{2} e^{in\\omega t} + \\frac{a_n + ib_n}{2} e^{-in\\omega t}], \\, 继续替换a_n,b_n 中的三角函数 \\\\ \\\\\n",
    "     & = \\frac{a_0}{2} + \\sum\\limits_{n=1}^{\\infty}[ \\frac{a_n - b_n}{2} e^{in\\omega t} + \\frac{a_n + ib_n}{2} e^{-in\\omega t}], \\, 继续替换a_n,b_n 中的三角函数 \\\\ \\\\\n",
    "     & = \\frac{a_0}{2} + \\sum\\limits_{n=1}^{\\infty}[\\, \\underbrace{\\frac{1}{T} \\int_{t_0}^{t_0+T} f(t) e^{-in\\omega t}\\,dt}\\, e^{in\\omega t} + \\underbrace{\\frac{1}{T} \\int_{t_0}^{t_0+T} f(t) e^{in\\omega t}\\,dt}\\,  e^{-in\\omega t}\\,] \\\\ \\\\\n",
    "     & = \\frac{1}{T}\\{ \\int_{t_0}^{t_0+T} f(t) e^{i0\\omega t}\\,dt\\,  e^{-in\\omega t}\\, + \\sum\\limits_{n=1}^{\\infty}[\\, \\int_{t_0}^{t_0+T} f(t) e^{-in\\omega t}\\,dt\\, e^{in\\omega t} ] \\, + \\sum\\limits_{n=-\\infty}^{-1}[\\, \\int_{t_0}^{t_0+T} f(t) e^{-in\\omega t}\\,dt\\, e^{in\\omega t} ] \\, \\} \\\\ \\\\\n",
    "     & = \\frac{1}{T} \\sum\\limits_{n=-\\infty}^{\\infty}[\\, \\int_{t_0}^{t_0+T} f(t) e^{-in\\omega t}\\,dt\\, e^{in\\omega t} ] \\\\ \\\\\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 黎曼积分  \n",
    "由欧拉公式，可以将$(1)$式写成如下形式：  \n",
    "$$f(t) = \\frac{1}{T} \\sum\\limits_{n=-\\infty}^{\\infty}[\\, \\int_{t_0}^{t_0+T} f(t) e^{-in\\omega t}\\,dt\\, e^{in\\omega t} ] \\tag{2}$$  \n",
    "如果我们要求一个函数$h(x)$在$x\\in[-\\infty,+\\infty]$上的积分$H(x)$，则有：  \n",
    "$$H(x) = \\lim\\limits_{\\Delta x \\rightarrow 0} \\sum\\limits_{x=-\\infty}^{\\infty}(h(x+\\Delta x) \\cdot \\Delta x) =  \\lim\\limits_{\\Delta x \\rightarrow 0} \\sum\\limits_{x=-\\infty}^{\\infty}(h(x) \\cdot \\Delta x) = \\int_{-\\infty}^{\\infty}f(x)\\,dx$$  \n",
    "令：  \n",
    "$$\\left\\{\\begin{array}{l}\n",
    " \\omega = \\frac{2\\pi}{T} \\\\\\\\\n",
    "\\omega_n = n\\omega \\\\\\\\\n",
    "\\mathcal{F}(\\omega_n) = \\int_{t_0}^{t_0+T}f(t)e^{-in\\omega t}\\,dt \\\\\\\\\n",
    "F(\\omega_n) = \\mathcal{F}(\\omega_n)e^{i\\omega_n t}\n",
    "\\end{array}\\right.$$\n",
    "$F(\\omega_n) = \\mathcal{F}(\\omega_n)e^{i\\omega_n t}$是一个关于$\\omega_n$的函数也是关于$\\omega$的函数，把$\\omega$看作$\\Delta \\omega_n$令其趋近于无穷小，则有：  \n",
    "$$\\begin{array}{ll}\n",
    "f(t) & = \\frac{T}{2\\pi} \\frac{1}{T} \\lim\\limits_{\\omega \\rightarrow 0} \\sum\\limits_{n=-\\infty}^{\\infty}[\\, F(\\omega_n) ] \\omega \\\\\\\\\n",
    "     & = \\frac{1}{2\\pi}  \\lim\\limits_{\\Delta\\omega_n \\rightarrow 0} \\sum\\limits_{n=-\\infty}^{\\infty}[\\, F(\\omega_n) ] \\Delta\\omega_n \\\\\\\\\n",
    "     & = \\frac{1}{2\\pi} \\int_{-\\infty}^{\\infty}\\mathcal{F}(\\omega_n)e^{i\\omega_n t}\\,d\\omega_n\n",
    "\\end{array}$$  \n",
    "\n",
    "其中：  \n",
    "$$\\left\\{\\begin{array}{llll}\n",
    "\\mathcal{F}(\\omega_n) & = & \\int_{t_0}^{t_0+T}f(t)e^{-in\\omega t}\\,dt && 傅里叶正变换 \\\\\\\\\n",
    "f(t) & = & \\frac{1}{2\\pi} \\int_{-\\infty}^{\\infty}\\mathcal{F}(\\omega_n)e^{i\\omega_n t}\\,d\\omega_n && 傅里叶逆变换 \n",
    "\\end{array}\\right. \\tag{3}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 拉普拉斯变换  \n",
    "拉普拉斯变换是建立在傅里叶变换的基础上的，为了解决函数不收敛的问题。强制加一个指数函数让大部分函数都能收敛。因为常用在信号处理、控制系统，其正逆变换公式如下：  \n",
    "\n",
    "$$\\left\\{\\begin{array}{llll}\n",
    "\\mathcal{L}(s) & = & \\int_{0}^{\\infty}f(t)e^{-s t}\\,dt && 拉普拉斯正变换 \\\\\\\\\n",
    "f(t) & = & \\frac{1}{2\\pi i} \\lim\\limits_{T \\rightarrow \\infty} \\int_{-\\gamma+iT}^{\\gamma-iT}\\mathcal{F}(s)e^{st}\\,ds && 拉普拉斯逆变换 \n",
    "\\end{array}\\right. \\tag{3}$$\n",
    "其中$\\gamma$是一个使$F(s)$的积分路径在收敛域内的实数。实际上使用会查表就够了！  \n",
    "\n",
    "#### 传递函数转化为微分方程    \n",
    "可以用微分算子$p=\\frac{d}{dt}$来替换方程中的$s$，并用时域函数$c(t),r(t)$替换复域函数$C(s),R(s)$，就能完成传递函数与微分方程间的转换。  \n",
    "$$\\begin{array}{l}\n",
    "G(s) = \\frac{C(s)}{R(s)} = \\frac{5}{s^2+3s+2} \\\\\\\\\n",
    "(s^2+3s+2)c(t)=5r(t) \\\\\\\\\n",
    "(p^2+3p+2)c(t)=5r(t) \\\\\\\\\n",
    "\\frac{d^2c(t)}{dt^2}+3\\frac{dc(t)}{dt}+2c(t)=5r(t)\n",
    "\\end{array}$$   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Python 代码  \n",
    "在Python 中直接求拉普拉斯的逆变换有时候会比较复杂:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 11 \\left(- 2 e^{2 t} + 5 e^{t} - 3\\right) e^{- 6 t} \\theta\\left(t\\right) + 3 \\left(e^{2 t} - 2 e^{t} + 1\\right) e^{- 6 t} \\theta\\left(t\\right) + 6 \\left(8 e^{2 t} - 25 e^{t} + 18\\right) e^{- 6 t} \\theta\\left(t\\right) + \\mathcal{L}^{-1}_{s}\\left[\\frac{s^{3}}{s^{3} + 15 s^{2} + 74 s + 120}\\right]\\left(t\\right)$"
      ],
      "text/plain": [
       "11*(-2*exp(2*t) + 5*exp(t) - 3)*exp(-6*t)*Heaviside(t) + 3*(exp(2*t) - 2*exp(t) + 1)*exp(-6*t)*Heaviside(t) + 6*(8*exp(2*t) - 25*exp(t) + 18)*exp(-6*t)*Heaviside(t) + InverseLaplaceTransform(s**3/(s**3 + 15*s**2 + 74*s + 120), s, t, _None)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *  \n",
    "from sympy.abc import a, s, t, omega\n",
    "G = ((s + 1)*(s + 2)* (s + 3))/((s + 4)*(s + 5)*(s + 6))\n",
    "inverse_laplace_transform(G,s,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 - 30/(s + 6) + 24/(s + 5) - 3/(s + 4)\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\delta\\left(t\\right) - 3 e^{- 4 t} \\theta\\left(t\\right) + 24 e^{- 5 t} \\theta\\left(t\\right) - 30 e^{- 6 t} \\theta\\left(t\\right)$"
      ],
      "text/plain": [
       "DiracDelta(t) - 3*exp(-4*t)*Heaviside(t) + 24*exp(-5*t)*Heaviside(t) - 30*exp(-6*t)*Heaviside(t)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "print(G.apart(s))  # 然后再查表或者分部逆变换\n",
    "inverse_laplace_transform(G.apart(s),s,t)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sympy 也不是万能的，很多时候简单的式子它也不一定能求得出来。还是记得展开为部分分式然后查表比较可靠！   \n",
    "或者利用`residue()`函数展开：\n",
    "$$G(s)=\\frac{0.1s^2+0.2s+1}{s^3+s^2+s}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.  +0.j        , -0.45+0.20207259j, -0.45-0.20207259j]),\n",
       " array([ 0. +0.j       , -0.5+0.8660254j, -0.5-0.8660254j]),\n",
       " array([], dtype=float64))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# s域部分分式展开 \n",
    "# r 留数（分子）\n",
    "# p 极点（分母）\n",
    "# k 常数项\n",
    "\n",
    "import scipy.signal as sig\n",
    "num=[0.1,0.2,1]\n",
    "den=[1,1,1,0]\n",
    "r,p,k=sig.residue(num,den)\n",
    "r,p,k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "结果为：  \n",
    "$$G(s)=\\frac{1}{s} + \\frac{-0.45+0.20207259j}{s-(-0.5+0.8660254j)} + \\frac{-0.45-0.20207259j}{s-(-0.5-0.8660254j)} + 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考  \n",
    "1. <https://github.com/alchemyst/Dynamics-and-Control/blob/master/1_Dynamics/3_Linear_systems/Laplace%20transforms.ipynb>  \n",
    "2. <https://zhuanlan.zhihu.com/p/41875010>  \n",
    "3. 《自动控制理论》  \n",
    "4. <https://tutorial.math.lamar.edu/pdf/Laplace_Table.pdf>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
